Two New Mitogenomes of Bibionidae and Their Comparison within the Infraorder Bibionomorpha (Diptera)

Despite the worldwide distribution and rich diversity of the infraorder Bibionomorpha in Diptera, the characteristics of mitochondrial genomes (mitogenomes) are still little-known, and the phylogenetics and evolution of the infraorder remains controversial. In the present study, we report complete and annotated mitogenome sequences of Penthetria simplioipes and Plecia hardyi representing Bibionidae. This is the first report of the complete mitogenomes for the superfamily Bibionoidea. There are 37 genes in each of the complete mitogenomes of all 20 studied species from eight families of four superfamilies within infraorder Bibionomorpha. The Ka/Ks analysis suggests that all 13 PCGs have undergone purifying selection. The gene rearrangement events exist in some families (Keroplatidae, Sciaridae, and Cecidomyiidae) but not in Mycetophilidae in Sciaroidea and also in Scatopsoidea, Anisopodoidea, and Bibionoidea, which suggests that these rearrangement events are derived in the late period in the evolution of the Bibionomorpha. The phylogenetic analysis suggests the phylogenetic relationships of Scatopsoidea + (Anisopodoidea + (Bibionoidea + Sciaroidea)) in Bibionomorpha. The divergence time analysis suggests that Bibionomorpha originated in the Triassic, Scatopsoidea and Anisopodoidea in the late Triassic, Bibionoidea in the Jurassic, and Sciaroidea in the Jurassic to the Cretaceous. The work lays a base for the study of mitogenomes in Bibionomorpha but further work and broader taxon sampling are necessary for a better understanding of the phylogenetics and evolution of the infraorder.


Introduction
The infraorder Bibionomorpha is one of the largest groups in Diptera with approximately 16,500 known species worldwide [1]. The insects in the infraorder are mostly tiny flies with a slender dark-colored body, long legs, and simple wings. Species in the infraorder are widely distributed across the world, and usually aggregate under moist shady areas, including leaves, grass clippings, moss, and manure. The feeding habits of their larvae are diverse, including detritophagy, saprophagy, predatory, mycophagy, and phytophagy [2]. Some larvae of bibionids and sciarids are found mostly in decaying organic materials such as forest litter, manure, and humus-rich soils, and play an important role as decomposers. Nevertheless, some species are known as economically significant pests, which cause damage to the roots of cultivated plants. The larvae of Bradysia cause serious damage to onions, carrots, and edible mushrooms [3], and their adults serve as vectors for the pathogen Fusarium foetens on the ornamental plant begonias, causing serious economic loss [4,5].
The classification of the Bibionomorpha has been studied for almost 200 years, and its monophyly has been confirmed by morphological and molecular studies [6][7][8]. The earliest comprehensive classification of the infraorder included 10 families primarily based

Sample Collecting and Mitogenome Sequencing
Samples of two species representing Bibionidae, P. simplioipes and P. hardyi, were collected from the Wuyi Mountains of Fujian Province, China in October 2018, and were preserved in individual vials in silica. After morphological identification using traditional method [32], they were stored in 100% alcohol, and housed at −20 • C until the DNA extraction. Total DNA was extracted from thoracic muscle tissues using the DNeasy tissue kit (Qiagen, Hilden, Germany) following the manufacturer's instructions. Concentration of extracted genomic DNA was determined by Qubit 2.0 (Invitrogen, Shanghai, China). The genomic DNA library was constructed with 350 bp small fragment, and sequenced using Illumina HiSeq X sequencing technology in HuiTong Company at Shenzhen, China. The adapters and unpaired, short, and low-quality reads were removed from the raw dataset during the quality control by FastQC v 0.11.3 [33].

Mitogenome Assembly and Annotation
The clean reads of mitogenome were separated and assembled using de novo assembler SPAdes v 3.11.0 with default parameters using the known dipteran mitogenome sequences of Clephydroneura sp. and Drosophila yakuba as references [34]. The contigs of mitogenome were extracted and assembled into mitogenomes through searching against the reference sequences using PRICE (paired-read iterative contig extension) by NOVO-Plasty version 2.6.2 [35]. The position and sequences of PCGs, tRNAs, rRNAs, and CR were originally annotated using the MITOS web server [36], and then determined in comparison with published homologous mitogenome sequences in phylogeny-close species using MEGAX [37]. The secondary structures of tRNAs were predicted by the MITOS, ARWEN 1.2, and tRNAscan-SE web Server [38]. The annotation of the mitogenomes was corrected manually using the Geneious v 4.8.5 [39], and final mitogenomes were submitted to the GenBank database.

Mitogenome Characteristics Analysis
The complete mitogenome for two newly sequenced Bibionoidea species were cyclized using the CGView online server with default parameters [40] (Table 1). Base composition, codon usage, relative synonymous codon usage (RSCU), and amino acid content of 20 species of mitogenome were computed with MEGA v.7.0.26 [41]. Nucleotide compositional bias was calculated using the following formulae: AT-skew = (A − T)/(A + T) and GC-skew = (G − C)/(G + C) by MEGA v.7.0.26, and three-dimensional scatterplots of AT-Skew, GC-Skew and AT% were drawn using Origin Pro v.9.0 [42]. Selection pressure of the 13 PCGs was analyzed by calculating Ka (non-synonymous) and Ks (synonymous substitution) values with DnaSP v6.12.03 [43], and visualized using Origin Pro v.9.0. The gene rearrangement events of tRNAs and PCGs were identified in comparison with the gene order of the ancestor Clephydroneura sp. using the Geneious v 4.8.5.

Phylogenetic Analysis
Phylogenetic analysis was conducted with both ML and BI using PhyloSuite v.1.2.2 [44]. Two datasets of mitogenomes were used for phylogenetic inference, and two species from the superfamily Tipuloidea, Symplecta hybrida and Tipula cockerelliana were employed as outgroup ( Table 2). we constructed two nucleotide datasets: (1) PCG123, all codon positions, and (2) PCG123 + R, all codon positions + two rRNA genes. The nucleotide sequences of 13 PCGs were aligned by codon-based multiple alignments using the L-INS-i algorithm and the rRNAs were aligned using the Q-INS-i strategy in MAFFT [45]. The concatenation of aligned sequences was performed using SequenceMatrix [46]. The selection of bestfit partitioning schemes and substitution models for each dataset were calculated using PartitionFinder 2 with the following settings: branch lengths as linked, and model election as AICc with the greedy algorithm [47]. Partitioning schemes and models are listed in Table S1. ML analysis was performed using IQ tree ver.1.6.8 under edge-linked partition model with 1000 bootstrapped replicates (BP) [48,49]. BI analysis was conducted using MrBayes ver. 3.26 [50] with two MCMC runs, each with four chains (three heated and one cold) run for 10,000,000 generations. Posterior probabilities (PPs) were computed after a 25% burn rate and each set was sampled every 1000 generations. An average deviation of the split frequency of less than 0.01 indicates that the runs had converged. The phylogenetic tree was visualized using FigTreev.1.4.4 and iTOL online tool [51,52].

Divergence Time Estimation
The divergence time was estimated using BEAST v.1.8.4 [53]. A γ distribution (GTR + I + G) nucleotide substitution model was selected by PartitionFinder 2.0 using AIC, and the speciation Yule model were selected as the tree priors with the uncorrelated lognormal relaxed molecular clock model. Two independent Markov chain Monte Carlo (MCMC) runs, each with a chain length of 100,000,000 generations with sampling every 1000 generations and a first 25% burn-in, were performed to estimate the divergence time. The fossil record of the species Eoditomyia primitiva 180 million years ago (Mya) for Sciaroidea was selected for calibration using the fossil calibration database and reported research searching in this study [54].

Mitogenome Characteristics of Bibionomorpha
The complete mitogenomes of P. simplioipes (GenBank no.: MT511121) and P. hardyi (GenBank: MT511122) are both of circular, closed, and double-stranded structures, with full lengths of 15,349 and 15,456 bp, respectively ( Figure 1). The complete mitogenomes of 20 species contained 37 genes (including 13 PCGs, 22 tRNA genes, and 2 rRNA genes) and CR. There are 22 genes (9 PCGs and 13 tRNAs) located on the majority coding strand (J-strand), while the other 15 genes (4 PCGs, 9 tRNAs, and 2 rRNAs) are on the minority strand (N-strand). The length of 20 mitogenome sequences ranges from 14,503 bp (Rhopalomyia pomum) to 16,950 bp (Acnemia nitidicollis), and the length variation mainly results from the CR, intergenic overlap, and spacers. The nucleotide base compositions and the three-dimensional scatter plot of the AT content, AT-skew, and GC-skew of 20 mitogenomes in the infraorder are shown in Table S2 and Figure 2. They all display obvious AT bias with A + T content ranging from 75.8% (Allodia sp.) to 85.7% (Orseolia oryzae). AT-skew values range from −0.034 (Arachnocampa flava) to 0.105 (Orseolia oryzae), and GC-skew from −0.326 (Pnyxia scabiei) to −0.102 (Epicypta sp.). Most of the species of the Bibionidae, Cecidomyiidae, and Sciaridae have similar AT content and AT/GC-skew, which are closely distributed in the three-dimensional scatter plot, whereas the species of the Mycetophilidae is widely distributed in the plot for AT content, AT-skew, and GCskew. The total PCGs nucleotide length ranges from 10,794 bp (Orseolia oryzae) to 11,237 bp (Sylvicola fenestralis). Most of the PCGs initiate with the typical start codon ATN and TTG, whereas the special start codons TCA, TTG, and TCG are found for cox1; TTG for nad1; and GTG for nad5. The most frequently used stop codons are TAA and T, followed by the stop codons TAG and TA. The RSCU values and amino acid usage frequency of the 20 species of mitogenomes in the infraorder are presented in Table S3 and Figure 3. UUA is the most frequently used codon, followed by UCU, CGA, and GCU. Leu has the highest usage percentage for all mitogenomes investigated with an average of 16.46%, and Cys has the lowest percentage (0.66%). The usage percentages of amino acids reveal no obvious difference among different families. To evaluate the selective pressures of 13 PCGs of the infraorder, pairwise analyses of the non-synonymous (Ka) and synonymous (Ks) substitution ratio (Ka/Ks) are shown in Figure 4. The Ka/Ks ratios for all 13 genes are all less than 1, ranging from 0.03 in cox1 to 0.98 in atp8 in the following order: cox1, cytb, cox3, cox2, atp6, nad5, nad3, nad1, nad4, nad2, nad4L, nad6, and atp8. These results imply that all of these 13 PCGs experienced purifying selection. Most tRNAs can be folded into the typical clover-leaf structure containing four stems and loops except trnS2, in which the dihydrouracil (DHU) arm is absent. The lack of a DHU arm in this infraorder has been commonly observed across all of these 20 mitogenomes. rrnL is located between trnL1 and trnV, and rrnS between trnV and CR. There are 12 mismatched base pairs (G-U) to be found in P. simplioipes tRNAs, and 13 mismatched base pairs (G-U) in P. hardyi ( Figure S1).

Gene Rearrangement Events of Bibionomorpha
Compared with the putative ancestral mitogenome (e.g., Drosophila yakuba), there are no gene rearrangements to be found in the superfamilies Scatopsoidea, Anisopodoidea, and Bibionoidea, and the family Mycetophilidae in Sciaroidea ( Figure 5); however, the

Gene Rearrangement Events of Bibionomorpha
Compared with the putative ancestral mitogenome (e.g., Drosophila yakuba), there are no gene rearrangements to be found in the superfamilies Scatopsoidea, Anisopodoidea, and Bibionoidea, and the family Mycetophilidae in Sciaroidea ( Figure 5); however, the

Gene Rearrangement Events of Bibionomorpha
Compared with the putative ancestral mitogenome (e.g., Drosophila yakuba), there are no gene rearrangements to be found in the superfamilies Scatopsoidea, Anisopodoidea, and Bibionoidea, and the family Mycetophilidae in Sciaroidea ( Figure 5); however, the arrangement events exist in some families in the Sciaroidea. In Keroplatidae, the trnE has been inverted in Arachnocampa flava. In Cecidomyiidae, the arrangement was observed in the genome upstream of nad5 (13 tRNAs and 5 PCGs), and downstream of nad5 with only one rearrangement (trnT and trnP). In Rhopalomyia pomum, three tRNAs (trnE, trnT, and trnP) have been inverted, trnN has been transposed from the position (trnR-trnS1) to a position (trnG-nad3), and trnI has been remote-inverted from the position (rrnS-trnQ) to a position (trnF-trnS1).
Genes 2023, 14, x FOR PEER REVIEW 9 of 18 arrangement events exist in some families in the Sciaroidea. In Keroplatidae, the trnE has been inverted in Arachnocampa flava. In Cecidomyiidae, the arrangement was observed in the genome upstream of nad5 (13 tRNAs and 5 PCGs), and downstream of nad5 with only one rearrangement (trnT and trnP). In Rhopalomyia pomum, three tRNAs (trnE, trnT, and trnP) have been inverted, trnN has been transposed from the position (trnR-trnS1) to a position (trnG-nad3), and trnI has been remote-inverted from the position (rrnS-trnQ) to a position (trnF-trnS1). In Mayetiola destructor, three tRNAs (trnE, trnT, and trnP) have been inverted, two tRNAs (trnN and trnY) have been transposed from typical position to the position (trnG-nad3) and (trnS1-trnE), respectively, and trnI has been remote-inverted (i.e., translocation + inversion) from a typical position to the position trnR-trnS1. In Orseolia oryzae, ten tRNAs have been remote-inverted, and five tRNAs and five PCGs have been transposed. In Sciaridae, the trnC has been remote-inverted from the position (trnW-trnY) to a position (trnQ-trnM), and the trnA and trnR were switched in Trichosia lengersdorfi. The trnC has been remote-inverted from the position (trnW-trnY) to the position (trnQ-trnM), and the trnY and the trnL2 have been remote-inverted, both of them from between trnC and cox2 to between rrnS and nad2 in Pseudolycoriella sp.

Phylogenetic Relationships and Divergence Time
Four phylogenetic trees were generated based on datasets PCG123 + R and PCG123 of 24 mitogenomes in Bibionomorpha. Two same topologies of trees from BI analysis ( Figure 6) are slightly different from the two same trees from ML (Figure 7) in the positions Arachnocampa flava in Keroplatidae and Pnyxia scabiei in Sciaridae. The former is grouped with Cecidomyiidae in BI, but at the base of Cecidomyiidae + Sciaridae in ML, and the latter with Trichosia lengersdorfi in BI, but at the base of Sciaridae in ML.    The superfamily Scatopsoidea is located at the base of the infraorder Bibionomorpha, and the superfamily Anisopodoidea at the base of the superfamilies Bibionoidea + Sciaroidea. The Bibionoidea + Sciaroidea appears monophyletic with PP = 1 and BP = 100 (Figures 6 and 7), and both of them are sister groups. In Bibionoidea, there are only three species to be investigated, and monophyly of the family Bibionidae cannot be determined. In Sciaroidea, there are four families Mycetophilidae, Keroplatidae, Cecidomyiidae, and Sciaridae to be investigated. The family Mycetophilidae seems monophyletic (PP = 1 and BP = 100), and is a sister group with Keroplatidae + Cecidomyiidae + Sciaridae (PP = 1 and BP = 100). The Sciaridae are monophyletic with PP = 1 and BP = 100. The position of Keroplatidae cannot be determined due to the position difference from both BI and ML analyses.
The PCG + R dataset was used to estimate divergence time because PCG + R and PCG resulted in the same topologies of phylogenetic trees using BI, but PCG + R had higher node support values than PCG in the initial phylogenetic assessment. We choose Sciaroidea as the fossil calibration point, and it was reported to originate at 180 Mya by Grimaldi and Engel. The Bibionomorpha was inferred to originate at 218 Mya, at least, in the Triassic (Figure 8). The superfamily Scatopsoidea is the earliest-appearing group in the investigation in the infraorder, and is estimated to originate at 218 Mya in the late Triassic. Subsequently, the Anisopodoidea appeared at 208 Mya in the late Triassic, and Bibionoidea +

Mitogenome Organization and Characteristics
The present study sequenced and annotated mitogenomes of two species in Bibi onidae, which is the first report of complete mitogenomes for the superfamily Bibi onoidea. There are 37 genes (13 PCGs, 22 tRNAs, and 2 rRNAs) and a CR, with high AT bias in nucleotide composition; the UAA, CGA, and GGA are three most used codons and the Leu, Phe, and Ile are the most used amino acid in 20 species of Bibionomorpha mitogenomes investigated in this study. All tRNAs have a complete secondary structure except for trnS2, which lacks the DHU arm. All of these characteristics are identical to the reported mitogenomes of other dipterans [27,55]. The Ka/Ks ratios in 13 PCGs range from the 0.03 in cox1 to the 0.98 in atp8, which implies that these PCGs have undergone purify ing selection. This is consistent with earlier studies in Culicidae [55,56], and in Mycetophi lidae [57].

Mitogenome Organization and Characteristics
The present study sequenced and annotated mitogenomes of two species in Bibionidae, which is the first report of complete mitogenomes for the superfamily Bibionoidea. There are 37 genes (13 PCGs, 22 tRNAs, and 2 rRNAs) and a CR, with high AT bias in nucleotide composition; the UAA, CGA, and GGA are three most used codons, and the Leu, Phe, and Ile are the most used amino acid in 20 species of Bibionomorpha mitogenomes investigated in this study. All tRNAs have a complete secondary structure except for trnS2, which lacks the DHU arm. All of these characteristics are identical to the reported mitogenomes of other dipterans [27,55]. The Ka/Ks ratios in 13 PCGs range from the 0.03 in cox1 to the 0.98 in atp8, which implies that these PCGs have undergone purifying selection. This is consistent with earlier studies in Culicidae [55,56], and in Mycetophilidae [57].

Gene Rearrangement of Mitogenomes
There are some rearrangement events in a small number of species in the superfamily Sciaroidea, but not in analyzed representatives of Scatopsoidea, Anisopodoidea, and Bibionoidea. The rearrangements exist in the families Keroplatidae, Sciaridae, and Cecidomyiidae but not in Mycetophilidae in the Sciaroidea. The unique inversion in Arachnocampa flava in Keroplatidae has been reported in other dipteran species, such as Anopheles quadrimaculatus and An. gambiae in Culicidae [58,59]. These three inversions, seven transpositions, and fourteen inverse transpositions (remote inversion) in three species in Cecidomyiidae have also been reported in Paracladura trichoptera in Trichoceridae in Diptera [11]. These two inverse transpositions and three transpositions in Trichosia lengersdorfi and Pseudolycoriella sp. in Sciaridae confirm earlier reports in Sciaridae [25].
Our phylogenetic study shows that the Sciaroidea is located at the top of four superfamilies investigated in this study, and Mycetophilidae and Keroplatidae + Cecidomyiidae + Sciaridae are sister groups in the Sciaroidea. These rearrangements exist only in the clade Keroplatidae + Cecidomyiidae + Sciaridae, which suggests that these rearrangement events appeared in the late period in the evolution of the Bibionomorpha. Their evolutionary mode and dynamics need to be explored with more mitogenome studies.

Phylogenetics and Evolution
The present study suggests the phylogenetic relationships of Scatopsoidea + (Anisopodoidea + (Bibionoidea + Sciaroidea)) in Bibionomorpha. Earlier molecular-based studies proposed Anisopodoidea to be thet most primitive group in Bibionomorpha [8,10], and another molecular-based study even suggested that the Anisopodoidea be included in one other infraorder Ptychopteromorpha [11]. Our result shows that the superfamily is derived later than Scatopsoidea, but only one species is investigated for both Anisopodoidea and Scatopsoidea in our study. Therefore, the phylogenetic position and monophyly of these two superfamilies are not yet settled. Although based on a very limited taxon sampling, this study supports the earlier research results for the monophyly of Sciaroidea [10,15]. The monophyly of Bibionidae in Bibionoidea is not yet determined due to only two species being included. The present study supports the monophyly of Mycetophilidae and Sciaridae in Sciaroidea, which was reported also based on mitogenomes [25,57]. However, further work and broader taxon sampling are necessary.
The present study suggests the origin of Bibionomorpha in the Triassic, which is consistent with an earlier study on Diptera evolution [60][61][62][63]. During the Triassic period, there were many autodigested fungi and decaying matter, which provided food for early lineages of the infraorder. The superfamilies Scatopsoidea and Anisopodoidea are suggested to be derived in the late Triassic and be the most primitive groups in the infraorder, which is consistent with the results of the previous study [9]. The Bibionoidea diverged in the Upper Jurassic (149 Mya), which is consistent with an earlier study [64]. Keroplatidae dates back to the Jurassic (170 Mya); the inferred dates are earlier than the inference by Blagoderov (105.3-99.7 Mya) [65]. The Mycetophilidae (150 Mya) and Cecidomyiidae + Sciaridae (155 Mya) in Sciaroidea are proposed to be derived in the Jurassic for the first time. The families Sciaridae (134 Mya) and Cecidomyiidae (101 Mya) in Sciaroidea originated in the Cretaceous, which is consistent with anearlier study [8]. In further research, more accurate estimates of divergence times are necessary with more precise fossil records for calibration and more complete sampling.

Conclusions
This is the first comprehensive study on the general characteristics of mitogenomes and mitogenomes-based phylogenetics in Bibionomorpha. Two species of mitogenomes are newly sequenced and annotated, and they represent the first report of complete mitogenomes for the family Bibionidae. A total of 20 mitogenomes in Bibionomorpha have been analyzed in the present study, and they are of the same general characteristics as other Diptera insects. Our results show that there are no gene rearrangements found in the superfamilies Scatopsoidea, Anisopodoidea, and Bibionoidea, and the family Mycetophilidae in Sciaroidea, but gene rearrangements exist in some families (Keroplatidae, Sciaridae, and Cecidomyiidae) in the Sciaroidea, which suggests that these rearrangement events are derived in the late period in the evolution of the Bibionomorpha. Our results support the phylogenetic relationships of Scatopsoidea + (Anisopodoidea + (Bibionoidea + Sciaroidea)) in Bibionomorpha, and the monophyly of the families Sciaridae and Mycetophilidae. In addition, we inferred that the Bibionomorpha originated in the Triassic, Scatopsoidea and Anisopodoidea in the late Triassic, Bibionoidea in the Jurassic, and Sciaroidea in the Jurassic to the Cretaceous. Our results demonstrate that the mitogenomes can be used to infer phylogenetic relationships in Bibionomorpha to help understand the history and evolutionary biology of the infraorder Bibionomorpha.

Conflicts of Interest:
The authors declare no conflict of interest.